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Abstract. This article presents a study of the sufficient accuracy of post- 
Newtonian and numerical relativity waveforms for the most demanding usage case: 
parameter estimation of strong sources in advanced gravitational wave detectors. 
For black hole binaries, these detectors require accurate waveform models which 
can be constructed by fusing an analytical post-Newtonian inspiral waveform with 
a numerical relativity merger-ringdown waveform. We perform a comprehensive 
analysis of errors that enter such "hybrid waveforms" . We find that the post- 
Newtonian waveform must be aligned with the numerical relativity waveform to 
exquisite accuracy, about 1/100 of a gravitational wave cycle. Phase errors in 
the inspiral phase of the numerical relativity simulation must be controlled to 
< O.lrad. (These numbers apply to moderately optimistic estimates about the 
number of GW sources; exceptionally strong signals require even smaller errors.) 
The dominant source of error arises from the inaccuracy of the investigated 
post-Newtonian Taylor-approximants. Using our error criterium, even at 3.5- 
th post-Newtonian order, hybridization has to be performed significantly before 
the start of the longest currently available numerical waveforms which cover 30 
gravitational wave cycles. The current investigation is limited to the equal-mass, 
zero-spin case and does not take into account calibration errors of the gravitational 
wave detectors. 
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1. Introduction 

Coalescing black-hole binaries are amongst the most promising sources for the current 
and future gravitational wave (GW) detectors such as LIGO, Virgo, LCGT [Hill 1311] 
and LISA [HEllT]. After several years of instrument upgrades, the next generation of 
ground-based interferometers should be operational within five years (~ 2015-2016). 
At advanced sensitivities, currently predicted event rates for stellar mass binary black 
holes (BBH) with total masses < IQQMq range from 0.4 to 1000 per year (with 20 
being the "realistic" number given in [S]) detectable up to distances of ^several Gpc. 
(Henceforth, for simplicity, we consider LIGO to be representative of all ground-based 
GW detectors.) In contrast, LISA will be sensitive to massive BBH systems with 
individual masses of 10^ — Mq up to redshifts of ^20. 

In order to detect GWs and derive the emitting sources' physical properties, 
accurate source modeling of predicted GWs is required in the form of a vast family of 
waveform templates. Analytical weak-field approximation methods such as the post- 
Newtonian (PN) expansion in general relativity accurately describe their inspiral prior 
to merger, whereas numerical relativity (NR) is used to model the merger of the two 
bodies in the strong-field regime. Both techniques have come to fruition over the 
last decade, e.g. the reviews (HI EHl Ell EH- Because of the computational cost of 
numerical simulations (^ 100, 000 CPU-hours for long inspirals of generic precessing 
systems), current NR simulations cover < 15 orbits. These NR waveforms are then 
fused together with a waveform covering the earlier stage of the inspiral, calculated 
using high order PN expansions. Such fused waveforms are referred to as hybrid 
waveforms. 

Hybrids play an important role in the construction of PN inspired phenomenolog- 
ical jl3| 114] waveforms, which are intended for event detection. In addition, hybrids 
may guide the construction of easy-to-implement waveforms that are being developed 
for measurement purposes. Within the Ninja project [151 116] . moreover, GW data 
analysts search for hybrid waveforms embedded in detector noise of ground-based de- 
tectors to investigate efficiency of GW data-analysis pipelines. Finally, the Numerical 
Relativity and Analytical Relativity collaboration (NRAR) [T7] aims to construct ana- 
lytical waveform models that span the entire parameter space of spins and mass ratios. 
Such models then allow us to compute waveforms with little computational expense 
for any choice of parameters. Some of the waveform models pursued in this collab- 
oration might involve hybrid waveforms. All these applications rest on sufficiently 
accurate hybrid waveforms, motivating our study. 

Besides hybrid waveforms, analytical waveform models based on the effective-one- 
body (BOB) formalism are common, e.g. |18[ [T5]. BOB waveforms are a resummed 
extension of PN approximants, e.g. [ID]. In this approach, one fits the BOB model 
directly against NR waveforms, without the intermediate step of constructing hybrid 
waveforms. Such fits, too, will be sensitive to errors in the NR waveforms. Conclusions 
about required accuracy of NR waveforms for hybrid waveforms may inform required 
accuracies for the NR waveforms that arc used when fitting BOB models. 

Both detection and measurement require hybrid waveforms to be of sufficient 
accuracy, to avoid a negative impact on the signal processing in noisy detector data. 
In conjunction with earlier work [H], Lindblom et al. [32] laid the foundations for 
accuracy standards for model waveforms. The formalism is applicable to all model 
waveforms (including PN only, NR only, and hybrids). In addition, though not 
considered in this work, Refs. [23J1M1115] extend the original model waveform analysis 
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of Ref. |22j by incorporating instrument calibration errors and deriving more refined 
accuracy requirements witliin the time domain. 

Based on Ref. [22] and complementary to recent investigations [26l [HI Ell [28] . 
we implement waveform accuracy requirements that preserve detection efficiency and 
parameter estimation performance. We place bounds on the norm of the modeling 
error Shmodci = ^modci — /loxact, where /loxact is assumed to be some exact waveform 
(a solution of Einstein's equations at infinite precision), and /imodoi is an approximate 
model waveform. Specifically, the error measure ||(5/i||/||/i|| translates to a normalized 
noise-weighted cross-correlation of the modeling error Shnwdei with itself, and is thus 
instrument-dependent, but distance- independent. 

Our work computes this quantitative accuracy measure in order to assess different 
sources of error arising from hybrid waveforms for non-spinning, equal-mass binary 
black holes. Rather than emphasizing the immediate aim of event detection, we take a 
long-term perspective, and consider the ultimate accuracy requirements necessary for 
data-analysis once GW detections have become routine. The emphasis lies, therefore, 
on the optimal extraction of source characteristics from the observed waveforms. Due 
to the very nature of their construction, hybrids suffer from sources of error such as 
the type and location of matching the PN and NR waveforms, and the PN and NR 
waveforms themselves. For instance, PN methods rely on an expansion in powers of 
v/c (where v is the BBH's typical orbital velocity and c is the speed of light jj, and 
result in an inspiral waveform with both PN amplitude and PN phase corrections. In 
addition, several different PN inspiral waveforms may be computed at a particular PN 
order because of differences when truncating Taylor expansions (commonly referred 
to as Taylor Tl, Taylor T3, Taylor F2, etc., see e.g. Ref. [SHj). 

An immediate question that arises when constructing hybrids is where matching 
between their PN and NR waveforms should occur. This is of practical importance 
because of the computational expense incurred in producing NR waveforms. Stitching 
PN and NR waveforms as late as possible minimizes computational cost. However, 
PN waveforms, derived using approximate perturbative methods in general relativity, 
become increasingly inaccurate as merger is approached. We investigate here where 
best to match PN and NR waveforms. Other possible sources of error, all defined 
below and explored here within a unified framework, include numerical resolution, 
the extraction radius of NR waveforms, the accuracy and approximant of the PN 
waveform, the matching width, and the effect of error within the matching procedure 
itself. We also examine quantitatively how detection and measurement accuracies 
vary using different ground-based instrument noise curves, all of which have different 
regions of sensitivity. 

This work thus provides a comprehensive quantitative measure of error for 
waveforms from several sources. The paper is organized as follows: Sec. 11.11 describes 
related earlier work and Sec. 11.21 introduces the notation used in this paper. Sec. [5] 
describes possible errors in hybrid waveforms and how to quantify them. Sec.[3]outlines 
our methodology: we begin with an explanation of the PN and NR waveforms in 
Sees. 13.11 and 13. 2[ of the hybridization process in Sec. 13. 3[ and of the computation of 
our error measures in Sec. 13.41 In Sec. HI we present a detailed study of the different 
sources of error in hybrid waveforms. Sec. [S] presents our conclusions. 

I The notation IPN corresponds to the formal ~ 1/c^ level in a post-Newtonian expansion with 
respect to the Newtonian acceleration and gravitational wave flux (where c is the speed of light). 
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1.1. Previous work 

Earlier works have examined similar aspects of detection and measurement errors of 
hybrid waveforms in the case of Initial and Advanced LIGO. Although the literature 
review presented here is far from complete, we chronologically outline four works 
which follow methodologies closely related to our study. As discussed in detail in 
Sec. m quantities in these works such as overlap, mismatch, faithfulness, effectualness 
and inaccuracy are related to our error measure |j(5/i||/||/i||. Throughout the rest of 
the paper, we will return to them to compare our results with their findings. 

First, Buonanno et al. [26j compute overlaps between different time- and 
frequency-domain PN waveforms and EOB waveforms for four non-spinning binaries 
with specific masses (the systems are (1.38,1.42) M©, (9.5,10.5) Mq, (10,1.4) Mq 
and (4.8,5.2) Mq). The overlaps are maximized either over time- and phase-shifts 
(faithfulness), or over time- and phase-shifts and masses (effectualness). Ref. [26] does 
not consider hybrid or NR waveforms per se. Instead, it uses EOB waveforms that were 
calibrated with NR simulations. In the case of Initial, Enhanced and Advanced LIGO, 
Ref. concludes that the majority of Taylor approximants (excluding TaylorT3 and 
TaylorEt) are sufficient for detection purposes for masses < 12M0 whilst advocating 
the use of EOB for masses > 12M0. 

Second, whilst constructing a new set of phenomenological waveforms for aligned 
and anti-aligned spinning BBHs, Santamaria et al. |14| compute the normalized 
distance squared (||<5/i||/||/i||)^ (a quantity defined as the inaccuracy, see Ref. [H]) 
for hybrids constructed by matching PN and NR waveforms in the frequency domain. 
This differs from our study, where matching proceeds only in the time domain. The 
PN waveforms comprise Taylor Tl, Taylor T4, or Taylor F2, and the NR waveform in 
the non-spinning case is an earlier version of the SpEC waveform used here. Santamaria 
et al. |14| determine the matching region itself by means of a least-squares fit between 
PN and NR waveform. This results in a matching region which extends very close 
to merger, Mu) G [0.063,0.126]. Quantifiable error sources are given by numerical 
resolution and by the differences between Taylor Tl, Taylor T4 or Taylor F2 hybrids. 
Ref. [13] concludes that the PN truncation error, represented by the difference between 
different PN approximants, provides the largest source of error. 

Third, and most similar to the approach followed here, Hannam et al. |27j aim 
to determine the necessary length of NR waveforms for the construction of hybrid 
waveforms. Ref. |27j considers both the mismatch (minimized over time and phase- 
shifts with or without mass-optimization), equivalent to one minus the overlap, and 
the inaccuracy, (||(5/i||/||/i||)^. As in Ref. [13], Ref. [37] constructs hybrid waveforms for 
aligned- and anti aligned- spinning and non-spinning equal- and unequal-mass BBHs. 
Focusing on the non-spinning equal-mass case, they use PN Taylor Tl and Taylor 
T4 approximants, and an earlier variant of the NR waveform used here. Ref. [57] 
computes quantitative error measures by comparing hybrids constructed from Taylor 
Tl and Taylor T4 with different matching frequencies between 0.045 and 0.09, using 
an Advanced LIGO noise curve with a low- frequency cut-off at 20IIz. Ref. [27] focuses 
on the near-term goal of event detection, and therefore places significantly weaker 
accuracy constraints on the hybrid waveforms. Nevertheless, Ref. [57] concludes that 
PN variants and PN order contribute the largest sources of error. 

Finally, Damour et al. [28] compute the effectualness and inaccuracy between PN, 
phenomenological and EOB waveforms for non-spinning unequal-mass systems. The 
PN waveform used is Taylor F2, a frequency domain approximant also employed in 
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Santaman'a et al. [H]. Similar to Ref. they do not consider NR waveforms 
and use an analytical Advanced LIGO noise curve [13]. Indeed, although similar 
in approach, Ref. [28] investigates different sources of errors in hybrid and EOBNR 
waveforms than the present work. Ref. 1281 finds that frequency-domain PN hybrids 
and phenomenological waveforms, unlike their EOBNR counterparts, fail to satisfy the 
minimal required accuracy standard for total mass M < 75Mq, even for moderately 
low signal-to-noise ratio. 

Shortly after the present article was submitted, a further study was performed 
by Boyle [30] . This reference studies length requirements of Taylor-PN approximants 
hybridized with EOB waveforms. 

Our work explores sources of errors in hybrid waveforms, aspects of which were 
investigated by these earlier works [5^1 dH [13 dH], within a uniform framework. 
In agreement with others, we find that the dominant source of error is due to the 
uncertainty between PN Taylor-approximants. Indeed, using our error criteria, even 
at the highest available PN accuracy, we show that NR waveforms must be longer 
than the currently available 30 GW cycles when constructing hybrids with Taylor PN 
approximants in the non-spinning equal mass BBH case. 

1.2. Notation 

Throughout the paper, we use the following notation: 

• hexact represents the exact waveform. 

• /ipN and /iNR denote the post-Newtonian and numerical relativity waveforms 
respectively. 

• /iH specifies a hybrid waveform. 

• p refers to the signal-to- noise ratio (SNR). 

• is our error measurement of the difference between a reference and trial 
hybrid. 

• O is the overlap between two waveforms. 

• is the mismatch between two waveforms, defined as 1 — C 

• pcff is the SNR divided by a detection safety factor e, and corrected for a network 
of detectors. Its inverse l/poff is the upper limit for our error calculations. 

2. Quantifying errors 

As mentioned in the introduction, hybrid waveforms are susceptible to a wide range of 
possible errors. We consider evolution of a binary black hole system with certain given 
intrinsic parameters like mass ratio, spins, and orbital elements. If one were able to 
solve Einstein's equations to infinite precision, one would obtain the exact waveform 
hoxact = hexactit) emitted by this system. For simplicity, we will only consider the 
dominant / = 2, m = 2 of the gravitational radiation. 

Numerical relativity computes an approximation to ft-oxact, which we shall denote 
/iNR- The NR waveform /inr niay differ from the exact waveform /loxact for a variety 
of reasons, among them: 

(i) Truncation error due to finite numerical resolution. As numerical resolution is 
increased, the solution will converge to the solution of the underlying continuous 
initial-boundary- value problem. 
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(ii) This continuous initial-boundary-value problem may differ from the realistic 
physical system, for instance, due to imperfect outer boundary conditions. 

(iii) Gravitational radiation extracted at a finite radius may differ from the real 
waveform at future null infinity, due to ambiguities in the definition of 
gravitational radiation at finite radius, and imperfect wave extraction algorithms. 

(iv) Initial data for the numerical simulation may not precisely correspond to the 
desired situation; for instance, the orbital eccentricity of the simulation might 
differ slightly from the desired value. 

This list is incomplete. We further emphasize that many possible sources of error for 
/iNR do not converge away with increased numerical resolution, such as, for example, 
the last three items in the list above. 

The most striking difference between /inr and /icxact is the length. The exact 
waveform is infinitely long, with an infinite number of gravitational wave cycles before 
merger. In contrast, the NR waveform has a finite number of cycles. Therefore, it 
is customary to attach a PN inspiral ft.pN at the beginning of /inr, resulting in a 
hybrid waveform hn of sufficient length to encompass the entire sensitivity band of 
the detector within a desired mass-range. 

Hybridization introduces several new potential sources of difference between ft-oxact 
and the hybrid waveform h^: 

(v) PN theory is in itself only an approximation to full general relativity; hence 
/iPN will differ from /loxact by some amount, which increases with decreasing 
separation of the binary. The accuracy of /ipN also depends on the PN order of 
the approximant. 

(vi) One must identify and match parameters like masses and spins in the PN 
formalism with those in the numerical simulatior(^. 

(vii) The matching between ft.pN and /inr is very sensitive to certain seemingly small 
errors in /ipN or /inr- 

To illustrate possible issues related to the matching, consider a particularly simple 
hybridization procedure. This procedure chooses a reference frequency a;o and time- 
shifts /iNR and /ipN such that both waveforms pass through this frequency at the 
same time. Here, we parameterize frequency by the PN-frequeney parameter [9] 
X = (M57)^/^, where M is the total mass of the binary and fl its orbital frequency. 
The sensitivity arises because the inspiral rate x decays very rapidly with decreased 
orbital frequency (i.e. increased separation), with i (x x^ during the inspiral at lowest 
PN order. Therefore, a small frequency-error dx in one of the two waveforms that are 
to be matched results in an error 

^i=-oc^ (1) 
X x° 

in the time-offset necessary to align the two waveforms with each other. The term 
x~^ causes the sensitivity of the matching between PN waveform and NR waveform. 
Matching typically occurs at frequencies x ~ 0.1. Whether or not St increases or 
decreases as one pushes the matching to lower frequency depends on how quickly Sx 
decreases with earlier matching. Let us consider the error due to the (unknown) 4PN 



§ This point was discussed in great detail in 1141 : here we simply choose masses and spins in the PN 
waveform to be equal to those in the NR simulation. 
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term contribution. As a function of timc-to-coalesccncc 0, the 4PN term contributes 
an amount 

Sx4PN OC Q-^^^ oc (2) 

to the frequency. The first proportionahty arises from the PN expansion of a; as a 
function of Q, (cf. [9]), whereas the second proportionahty arises from the Newtonian 
relation x oc 9~^/^. Substituting Eq. ^ into Eq. we see that the unknown 4PN 
coefficient contributes an error 5t independent of the matching frequency. 

As a second example, let us consider the orbital phase: the (unknown) 4PN term 
of the orbital phase, expressed as a function of frequency x is 

$4PN OC ^3/2 OC (3) 

The impact of unknown higher-order PN terms on the orbital phase, therefore, 
decays with smaller frequency x. However, this decrease is very slow: To reduce 
the effect of $4pn by a factor of 2, one has to terminate the PN waveform by a 
factor 2^/^ « 6.4 earlier in time, with a corresponding increase of the length of the 
NR simulation which takes the binary to merger. Because of the high computational 
cost of an NR simulation, increasing its length by large factors like 6.4 is hardly 
practica|i|. Furthermore, the numerical waveform /inr also contributes potential issues 
for matching at earlier times. Widely separated binaries are particularly difficult to 
simulate numerically, due to the small energy flux which drives the inspiral, and due 
to the long wavelength of the gravitational radiation. Moreover, the first portion of 
numerical simulations tends to be noisy due to junk radiation |34| I35[ I36| 137] arising 
from the initial conditions. 

We see, therefore, that there are a wide variety of effects which might cause 
to differ from the desired ft-exact • Whether such differences are relevant will depend on 
the application of the waveforms. We focus here on applications to GW detection and, 
in particular, measurement. The sensitivity of GW detectors is frequency dependent 
(see Fig. [T]), and therefore, the importance of errors will depend on their frequency 
dependence. Because BBH waveforms are invariant under a rescaling of mass, one 
PN-NR hybrid waveform /ih represents binaries of any total mass M for a given mass 
ratio. However, depending on the mass M, different portions of are in the most 
sensitive frequency range of the detector. Furthermore, waveforms are typically used 
in matched filtering, e.g. [SHIEI], a technique that allows for time- and phase-shifts 
of the gravitational waveform. Thus, those components of the difference /ih — ^exact 
that can be represented by a time- or phase-shift are not relevant. 

Let us assume that the exact signal /lexact and the computed hybrid waveform 
/iH are so close to each other that the Taylor-series expansions used implicitly in the 
Fisher-matrix formalism are valici^. We denote 

dh= hu- /icxact, (4) 

where Sh is the difference between the hybrid and exact waveform when optimally 
aligned in phase and time. To introduce some key concepts, we shall assume that we 
know 6h. Miller [5T] and Lindblom et al. [52] pointed out that the key quantity for 
assessing the importance of errors is the inner product = {6h, Sh), where 

II Implicit time-stepping techniques 1311 1321 l33l are expected to reduce this hmitation. 
^ A commonly made assumption in gravitational wave data-analysis. 
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Figure 1. LIGO noise curves used in our analysis. We predominantly use 
ZERO_DET_HIGH_P, as predicted for Advanced LIGO |40j. ZERO_DET_LQW_P is the 
Advanced LIGO sensitivity curve with a lower-powered laser. In Sec. 14.51 we 
compare to two earlier noise curves, Initial LIGO, and an earlier Advanced LIGO 
noisecurve used by Ajith et al. [13]. The wider bandwidth of ZERO_DET_HIGH_P 
places more stringent phase coherency on waveform templates. 



with g{f) and h{f) the Fourier transforms of two waveforms g(t) and h(t). Sn{f) 
denotes the (one-sided) power spectral density matrix, 

/oo 
dTe2--f-C„(r) , />0, (6) 
-oo 

where C„{t) is the noise correlation matrix for zero- mean, stationary noise. The 
noise spectra employed in this paper are shown in Fig. [Tj We will generally focus 
on the noise curve ZERO_DET_HIGH_P, the highest-sensitivity design goal for Advanced 
LIGO [10]. This configuration is sensitive over a wider bandwidth than the Initial 
LIGO, and earlier estimates for Advanced LIGO, used, for example, in Refs. [13] 
and [28]. Therefore, the FN inspiral and the NR late inspiral-merger-ringdown are 
simultaneously in band for a larger range of total binary masses than for previous 
noise curves. For the earliest matching frequency considered here {Mw = 0.038), the 
matching region will be in Advanced LIGO's sensitivity band (defined here as lOHz 
to lOOOHz) for total masses between about 1.2Mq and 125M0. 
When 

\\Sh\\ < 1, (7) 

a detector with noise-spectrum Sh{f) cannot experimentally distinguish between the 
waveforms ft-n and ft-oxact- Therefore, hn is a suitable substitute for /loxact for parameter 
estimation in this case. The inequality ([7]) is sufficient but not necessary. For instance, 
if Sh is orthogonal to the signal manifold, its impact on parameter estimation is 
diminished. 

Event detection genenerally places less restrictive demands on the accuracy of the 
waveforms hn. A sufficient criterion is given by Lindblom et al. [22) . 

\\Sh\\<V2^P, (8) 
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where the signal-to-noise ratio (SNR) p = |j/iexact||- The dimensionless number e 
depends on the reduction in SNR (and thus event rate) one is wilhng to tolerate 
during event detection. For typical template placement in Initial LIGO data-analysis, 
Lindblom et al. suggest Emax = 0.005. These numbers may have to be reconsidered 
for Advanced LIGO data analysis [H], based on the cost of using a template bank 
(which depends on the number of templates) relative to the cost of computing the 
templates in the first place. 

For non-spinning equal mass BBHs, the waveforms /loxact and hn depend trivially 
on two parameters: the luminosity distance and total redshifted mass. The luminosity 
distance to the source is introduced as an overall scaling of the amplitude. For 
BBHs (as considered here), the total redshifted mass comes into play only via trivial 
rescalings of amplitude and time (i.e. frequency). The dependence on distance and 
mass carries forward into \\6h\\, complicating evaluation of Eqs. ([7]) and (jS]). The 
distance dependence can be removed through division by the signal-to-noise-ratio, i.e. 
by computing 

\\Sh\\ 



P 

In terms of this quantity, the accuracy requirements becomcQ 



(9) 



\6h\\ { 1/Poff, parameter estimation 
1 1'* 1 1 1 \/2£max- event detection 



The square of the left-hand-side of Eq. (|10p is called inaccuracy functional by Damour 
et al. [28j . Because of our emphasis on numerical errors, we prefer to define this 
quantity without the square, so that ||(5/i||/||/i|| will be proportional to the numerical 
errors. This has the added benefit that the right-hand side is proportional to the 
inverse of the signal-to-noise ratio and not its inverse square. 

For the parameter estimation limit in Eq. (|10|) . we have replaced the single 
detector SNR p by an effective SNR pcflEl- This effective SNR allows us to incorporate 
a safety factor e < 1, as recently argued by Damour et al. [3S]. The safety factor 
reduces the impact of waveform errors 5h to a fraction e of the effects of the detector 
noise; Damour et al. [55] suggest 1/3 — 1/2, in light of the fact that it is easier 
to calculate more accurate waveform templates than to enhance the sensitivity of the 
GW detectors. Furthermore, p^s can absorb the impact of a network of detectors. 
In a network of detectors with Gaussian, independent noise, the overall network SNR 
is the root-square-sum of the detector SNRs. Thus, though GW detectors have non- 
Gaussian noise, one can still get a rough sense of the requirements for a coherent 
network of N individual detectors by incorporating a factor ^fN into peff, 

p,s^e-^^/Np. (11) 

Let us consider reasonable values for pcff- Our focus is on parameter estimation, 
i.e. we would like to obtain results that guarantee sufficiently accurate hybrids for 
most Advanced LIGO events. In realistic scenarios, Advanced LIGO may be seeing 
20 binary black holes per year [8]. Because LIGO is volume limited, most of these 

+ For event detection, it is irrelevant which template matches a certain physical signal hcxact, and 
reliable estimates of detection sensitivities require mismatch calculations optimizing over physical 
parameters like total mass or mass-ratio. We do not consider such optimizations in this paper, 
and therefore Eq. HlOj l gives only an approximate, conservative indication of requirements for event 
detection. 

* Lindblom et aJ.'s original bound can be recovered by sotting e = 1. 
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events will have an SNR close to the detection threshold, with the number of stronger 
events decaying like p~^. The strongest event in a year will have a SNR larger by a 
factor ~ 20"'^/'^ ~ 2.7 than the detection threshold, placing it around SNR of w 20. 
Including a safety factor e = 1/2, one arrives at poff = 40. Accounting for 3 detectors 
increases this bound to peS = 40-\/3 sa 70 and for 5 detectors to peS = 40-\/5 « 90. 
Using the optimistic event rate estimate of 1000/yr [S| instead of the realistic rate of 
20/yr, increases the largest expected SNR's by a further factor (1000/20)^/^ ~ 3.7 to 
PeB ~ 250. To cover this range of possibilities, we will indicate poff = 40 and pcs = 100 
in the plots below. In addition, the event detection limit of Eq. (fTO]). V2emax ~ 0.1, 
can be rewritten in terms of peff = 10. We also indicate this bound in our figures. 

The dependence of ||(5/i|| and ||5/i||/||/i|| on the total mass can be taken into account 
by plotting ||(5/i||/||/i|| as a function of total mass. Such a plot gives insight into 
suitability of hn for both event detection and parameter estimation; if ||(5ft,||/||/i|| is 
below V2emax ^ 0.1 for certain masses, criterion ^ is satisfied, and ft-n is suitable for 
event detection in that mass-range. The value of ||5/i||/||/i|| as a function of mass gives 
the (inverse) signal-to-noise ratio, up to which hn is suitable for parameter estimation. 
When a signal at a certain SNR pobs is observed, the plot allows one to verify whether 
||(5/i||/||/i|| < e/pohs- If this inequality is satisfied, event characterization can proceed 
without delay; otherwise, more accurate waveforms hn need to be computed for 
optimal event characterization. 

So far, we have assumed knowledge of the exact solution /lexact, which, in practice, 
one does not know. Instead, we will below compute differences Sh = g — h between 
two hybrid waveforms g and h, where we consider the "superior" of the two hybrid 
waveforms as a substitute for /icxact- Thus, we will usually consider 



Replacing the denominator of Eq. ^ by or by changes Eq. (fT^ 

only by terms of higher order in 6h and is therefore irrelevant. We will, therefore, 
consider Eq. ([T2|) as a measure of the error of the "inferior" of the two waveforms g, 
h. In doing so, it is important that one of the waveforms is clearly more accurate 
than the other, otherwise Eq. p2p would simply measure the difference between two 
poor-quality waveforms (and such a difference can be arbitrarily small, if the two 
poor-quality waveforms happen to be near each other). We will discuss this point in 
more detail below. 

The error measure ||(5/i||/||ft.|| can be expressed in terms of the overlap O and the 
mismatch M. = \ — O. The overlap between two waveforms g and h is defined as 



(.g,/i) _ {h + 5h\h) 



(13) 



^/{g,g){h,h) ^{h + Sh\h + Sh){h\h)' 



By Taylor expanding the denominator, we find: 




(14) 



and 




(15) 



Here, is the normalized projection of Sh onto h, 



(16) 
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Normalizing Eq. ((T5)) by instead of merely affects higher order terms. The 
quantity represents an overall resealing between h and g; ii g — {1 + for some 
constant 7, then =7. 

Eq. (|15p relates our error measure to the more widely used mismatch 

Ai. The major difference is that is also sensitive to an overall resealing of the 

waveform, represented by whereas Ai is not. Waveforms h are high-dimensional 
vectors either in the time-, or frequency-domain, and the quantity h\\ only measures 
one dimension of this high-dimensional vector space. If one assumes errors that have 
components in many dimensions, and consequently neglects , then Eq. ()15p simplifies 
to 

\\Sh\\/\\h\\ = V2M. (17) 

This relationship has also been noted in [22l [HI [271 112] • 

When we compute in subsequent sections, we always introduce a 

relative time- and phase-shift between g and h in order to minimize ||(5/i||/||/i|| and 
the mismatch A4 , but we do not modify any other parameters of the waveform (such 
as the total mass). Thus, our analysis tests faithfulness A more appropriate error 
limit for detection is the ejfectualness where the error is minimized over the intrinsic 
parameters of the system, that is, the total mass in an equal-mass, non-spinning 
binary. Thus, our upper error limit of p^f ~ 10 for detection is very conservative. 

Finally, we note that the analysis presented here is not the first one in this spirit. 
Mismatches or between various pairs of waveforms have been calculated 

before, for instance between purely numerical waveforms in the Samurai project [43], 
between analytical waveform models by Damour et al. [28], and also between hybrid 
waveforms |141 127j . Our focus is a comprehensive and unified study of many error 
sources that enter hybrid waveforms, concentrating on parameter estimation. 

3. Methodology 

3.1. Post-Newtonian waveforms 

The PN approximation in General Relativity relies critically on assumptions of 
weak field and of the source's internal slow motion. It is characterized by the PN 
parameter e ^ (w/c)^ ~ (M/r), where AI is the characteristic mass of the system 
and V the magnitude of the relative velocity. Construction of PN templates requires 
modeling both the local conservative dynamics of the system and the generation of its 
gravitational waves. The current state-of-the-art accuracy in both the gravitational 
wave generation and the description of the system's local dynamics is 3.5PN in the 
case of non-spinning quasi-circular comparable mass inspiralling compact binaries 
(see Blanchet [9] and references therein for a detailed review). In the case of 
spinning inspiralling compact binaries, the amplitude-complete GW templates have 
been computed to 1.5PN order [H]. Following standard notation, nPN describes the 
n-th PN order in the phase, and "PN approximant n/m" specifies the n-tli and m-th 
PN corrections in the phase and amplitude respectively. 

Standard formulae exist for PN gravitational polarization waveforms h^(t) and 
hx{t) of quasi-circular orbits, which depend on the source's location with respect to 
some observer. Usually, NR waveforms are decomposed into spherical harmonics and 
often only the dominant (2,2) mode is considered. This is the case with the set of 
NR simulations we use, and consequently, we will only need the (2,2) mode of the PN 
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2 /3 

where the PN invariant velocity parameter is x = {A-m{t)Y'-', n{t) is the orbital 
frequency and = Q,(t')dt' is the orbital phase. The parameter Inzp = 

^ — I7 + I In (^^) is a constant frequency scale which depends on the freely 
specifiable time parameter b. The latter characterizes the timescale at which the 
gravitational wave tail contributions enters the polarization waveforms. We choose 
6=1. 

Evaluation of (fTB]) requires expressions for x{t) and We will investigate 

the Taylor approximants TaylorTl, TaylorT3 and TaylorT4 (as defined in [46] and 
references therein) . They are all completed by solving the energy-balance equation: 

where the GW flux at infinity C is balanced by the change in the orbital binding energy 
E{t) of the binary. The TaylorT approximants differ from each other by how higher- 
order PN terms are truncated in E{t) and C when implementing In particular, 
TaylorTl and TaylorT4 arc solved by numerically integrating coupled sets of ordinary 
differential equations (see (35-36) and (45-46) in [IS] respectively and [3^1117]), whereas 
analytical expressions exist for TaylorT3 (see (43)- (46) in [3B] and [H]). Many previous 
studies have established that for equal-mass, non-spinning binaries the TaylorT4- 
approximant happens to agree very closely with the full numerical solution of Einstein's 
equations US] SHI SSI [50] ■ There is, however, no a priori reason why TaylorT4 should 
perform better than TaylorTl or TaylorT3. So that our results are as representative 
as possible, we therefore investigate several Taylor approximants (TaylorTl, T3, T4). 
By doing so, we achieve conservative estimates of how long NR simulations must be. 
We expect that our results are unaffected by accidental agreement of TaylorT4 for 
the specific configuration under study, so that our results will carry over to other 
parameter configurations where no exceptionally good PN approximant is known. 

We start the PN waveforms at an initial orbital frequency oi AdVl = 7.5 x 10""', 
resulting in a waveform of duration ^ 1.7 x l{f M containing ^ 6500 GW cycles to 
merger. The low starting frequency ensures an initial gravitational wave frequency 
of 10 Hz at M = 5M0, so that the constructed hybrid waveforms cover the entire 
frequency band of the Advanced LIGO detectors for M > 5AIq. 



3.2. Numerical Relativity Waveforms 

In this paper, we focus on equal-mass non-spinning BBHs. We use a waveform 
generated with the Spectral Einstein Code SpEC [51]. The initial data is identical 
to that used in |461 152j . resulting in an inspiral lasting about 16 orbits at an initial 
orbital eccentricity of ~ 5 x 10~^. The techniques for the evolution are a refinement 
of those presented in [IHJ [52] . Specifically, the gauge source functions are rolled off to 
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zero early in the evolution, so that the inspiral proceeds in the harmonic gauge. The 
constraint-damping parameters are those described in |53| . which reduce the impact 
of junk radiation. The merger was performed with a refinement of the techniques 
described in [54] , 



3. 3. Construction of hybrid waveform 

Construction of a hybrid waveform h}i{t) from a PN waveform hp^{t) and NR 
waveform h^Yiit) requires first that we align the waveforms with a relative time- and 
phase-shift, and then join them smoothly together. We begin by choosing a GW 
frequency interval of width Suj centered at frequency oj™, 

Sio 5uj 

in which we perform the matching between NR and PN. During the inspiral. the 
gravitational wave frequency uj{t) = d(f>{t)/dt is continuously increasing, so (PU)) 
translates into a time-interval imin < t < <max- Here, represents the phase of 
the (2,2) mode of the gravitational radiation. 

The PN waveform hp^{t;tc,^c) incorporates naturally the coalescence 
parameters tc and representing the time and phase freedom. We determine the 
parameters t'^ and $^ by minimizing the GW phase difference in the matching interval 

tc/$c = minfocy^ (0PN(t;tc,$c)-</'NR(t))'dt. (21) 

The hybrid waveform is then constructed in the form 

knit) = Ht)hpN{t; t'c, $^) + [1 - Fit)] h^nit), (22) 
where F(t) is a blending function defined as 




t <C ^min 

cos lll^l'^'^l ) , ^min < t < tmax (23) 

t ^ ^max 

Several other hybrid procedures have been proposed. Similar to our method, 
Hannam ct al. and Ajith ct al. |271I13| match in the time domain, although they employ 
different blending functions (e.g., Eq. (4.11) in [13]). Both Santamarfa et al. [14] and 
Damour et al. |28) match in the frequency domain using a fit, and invoke different 
blending functions. Earlier hybrids were matched at a single frequency j461 I50j . 
We have experimented with minimizing the phase and amplitude difference when 
determining time and phase offset. Specifically, we replaced the integrand of ^T\i by 
|/iPN — /iNRp using the complex (2,2) modes. We did not notice any significant change 
when using this different alignment scheme. 



3.4- Computation of \\Sh\\/\\h\\ 

We compute Fourier transforms of waveforms with standard Fast Fourier Transforms, 
rescaling the geometric units in which the hybrid waveforms are given to the desired 
physical mass. In order to reduce the Gibbs phenomenon, we multiply the time-domain 
waveform with a windowing function wit) before computing the Fourier transform. We 
tried two windowing functions, the Hann (or cosine squared) window |55| and the more 
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recent Planck-tapcr window function |56| . Our results did not differ noticeably when 
using the different windowing functions, and so the Planck-taper window |56] is used 
since it quickly goes to zero outside the window bounds: 

0, t<ti 
[ey(t) + i]~\ ti<t<t2 

W{t) = <( 1, t2<t<t3 (24) 

[e^(*) + l]"\ t3<t<t4 
0, ti < t, 

where 2/(0 = it2-ti)/it-ti) + {t2-ti)/{t-t2) and z{t) = -{t^-ti) / {t-ta) + {t^-ti) / (t-t^). 
The interval [ti , ^2] is chosen to cover the first five GW cycles of our hybrid waveform, 
whereas [^3,^4] spans the very late ringdown. 

Detector noise has a considerable effect on {h,g), since different parts of the 
waveform enter the detector band for different masses. When evaluating the noise 
weighted inner products in Eq. ([S]), we use the noise curves shown in Fig. [1] with the 
default being the Advanced LIGO ZERD_DET_HIGH_P curve. 

Whenever we calculate ||(5/i||/||/ij| between two hybrid waveforms, we time- and 
phase-shift one of the two hybrid waveforms to maximize the overlap for each mass. 
We can write the Fourier transforms of the time- and phase-shifted waveform as 
follows: 

hA4,,At = (he cos A(P + hs sin A(f>)e^'''f^* . (25) 

The extremization over the phase-shift A0 can be performed analytically. We define 
a = {hc,g) and b = {hs,g), where he is the Fourier transform of the real part of the 
first waveform, and hg is the Fourier transform of its imaginary part. Thus, 

A</> = tan^i Q^, (26) 

which is the phase shift for maximum overlap if [— sin A0 — cos Ac/)] < (otherwise 
TT is added to the phase). At is found by using the Matlab minimization function 
f minsearch. 



4. Results 

In the following error analysis, we focus primarily on the upper limits of ||(5/i||/||/ij 
for a single detector, corresponding to peS of 10 and 40 for detection and parameter 
estimation respectively. We do not discuss the upper error limit for multiple detectors, 
Peff ^ 100, since none of our hybrid waveforms meet this error criterion, with the 
possible exception of those hybrids matched very early with a TaylorT4 waveform. 

4-1- Effect of error in hybridization 

During the hybridization procedure, as outlined in Section fS.SI the FN waveform must 
be aligned with the NR waveform by choosing tc and <i>c- The time shift tc is difficult 
to determine, because a degeneracy between time and phase shift is only broken by the 
frequency evolution, which is very slow during the inspiral, cf. Sec. [5] On the other 
hand, minimization over the phase-shift $c is trivial, because $PN(i;tc,^c) — ^NR(i) 
in the integrand of Eq. (|2ip is linear in <f>c- Therefore, we assume that errors in 
the determination of $c are irrelevant, and focus on the time-offset by investigating 
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Figure 2. Effect of hybridization errors. We compare hybrid waveforms with 
alignment biased by ATbias to the correctly hybridized waveforms. (Matching 
frequency Muim = 0-050, highest-resolution NR waveform, and TaylorTS PN 
waveform) . 



how accurately tc must be determined. For this, we first construct a hybrid using 
the correct least-squares fit, resulting in the correct t'^. We then construct a second 
hybrid by fixing the time-ofi^set to -I- ATbias and by performing the minimization 
over $c only. We compute between these two hybrids. All least-squares 

fits are performed over a fi'equency range of width Suj ~ 0.1a;„i centered at frequency 
Mum = 0.05. 

Fig. [3] presents the results of this calculation. The error j|(5ft-||/||/i|| decreases 
linearly with ATbias; this linear dependence is one advantage of using ||(5/i||/||/i|| over 
other measures, such as its square or the mismatch. Fig.[2]shows that even small errors 
in hybridization have a significant effect on the magnitude of ||(5/i|j/||/i||. Time-shift 
errors larger than about 2M are unacceptable for parameter estimation at poff ^ 40. 
In order to preserve some safety margin, we therefore conclude that the time-shift 
must be determined to better than IM. (Strictly speaking, this estimate AT < IM 
applies only for matching at a frequency Muj„i ~ 0.05; this analysis would have to be 
repeated when matching earlier or later). This is a surprisingly tight bound; IM is 
only ~ 1/2000 of the time-to- merger, and is only about 1% of a gravitational wave 
cycle. 

Our primary goal is to utilize the bound on ATbias when choosing the width of 
the matching interval in Sec. 14.21 However, to place this requirement ATbias ^ 
into perspective, consider the TaylorT2 PN approximant, which expresses time-to- 
merger as a power-series in orbital frequency. The 3.5PN-order contribution to the 
coalescence time for an equal- mass binary is given by (see [9] and references therein); 

'-™^^^'<""'""-** 

where we have used Mfl = Mujm/'2. = 0.025. The contribution ia.spN is two orders of 
magnitude larger than the bound on ATbias- Whilst it is difficult to translate time- 
offsets from PN contributions, this large disagreement immediately indicates that PN 
theory, even at 3.5PN order, might not be sufficiently accurate. 
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Figure 3. Effect of tfio width of matcfiing window Sui/uJrn on the best-fit time- 
shift tc- For matching windows covering < 1 GW period {&u)/ujm % 0.02), tc is 
sensitive to small oscillations in the GW phase (thin lines) . For windows covering 
> 1 GW period, this sensitivity rapidly decreases (thick lines, also enlarged in the 
inset). 

We note that Hannam et al. [Uj also considered hybridization errors; this work 
compared three different hybridization schemes, and found that their difference is 
indistinguishable for certain SNRs. Our approach improves on this by explicitly giving 
conditions on the accuracy of <c as in Fig.[5J This measure of hybridization errors might 
seem extremely specific because we only consider the errors due to the determination 
of tc, however these error bounds apply to any time-domain hybridization procedure 
and to any effect that might give rise to an error in determination of tc, like for instance 
a small eccentricity in the NR simulation. 

4-2. Width of hybridization interval 

The PN time- and phase-parameters tc and $c that align the PN waveform with the 
NR waveform are determined by a least-squares fit [see Eq. (PT|) ]. The width of the 
matching interval [imin , ^max] influences the robustness of this fit for two reasons: 

(i) A small change Stc in tc is nearly degenerate with a small change S^c in ^c, with 
S^c ~ uj,nStc/2. (The factor 1/2 arises because w,„ measures the frequency of the 
GW I = 2, TO = 2 mode whereas $c is the orbital phase). For circular Newtonian 
orbits, the degeneracy is perfect. For the general relativistic case, the degeneracy 
is broken only by the evolution during the hybridization interval, i.e. by the 
change in frequency uj. A longer fitting interval encompasses more evolution, and 
so helps to break the degeneracy. 

(ii) A longer fitting interval helps to "average out" undesirable features in the NR 
waveform. Such features might be caused by numerical artefacts like junk 
radiation or a slight residual eccentricity. In particular, matching over an interval 
much smaller than the orbital period will be susceptible to eccentricity effects, as 
discussed in Boyle et al. [35]. 
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Figure 4. Effect of numerical errors. Shown is ||(5h||/||/i|| between hybrid 
waveforms constructed with NR waveforms of different numerical resolution 
(dashed lines) , and between hybrid waveforms of finite GW extraction radius /?gw 
vs. waveforms extrapolated to infinite extraction radius. (All hybrids matched at 
Muj = 0.042 to a TaylorTS PN waveform). 



A robust fit Eq. (j2ip must be independent of small changes to the fitting interval. To 
test this, we compute many independent least-squares fits, over slightly shifted fitting 
intervals, parametrized by the frequency in the middle of the fitting interval. 
Fig. [3] plots the resulting tc for several different widths 5uj of the matching interval. 

For small 5uj, this figure shows significant variations in t^. These variations arise 
because the fitting interval covers less than one GW period. Therefore, the fitting 
is sensitive to a small variation in the GW phase (of order 0.001 rad) present in the 
numerical data. As 5uj/u) > 0.05, this sensitivity rapidly decreases, as the fitting 
interval now encompasses multiple cycles, averaging over the small oscillations in the 
numerical waveform. The remaining overall trend (tc ~ 3660Af to tc ~ 3630Af ) is 
explained by the increasing error in the PN waveform at higher frequencies; the value 
~ 3650M arises because the time in the NR simulation is chosen such that t = at 
the beginning of the simulation, with merger around t ^ 4000Af. For Suj/oj = 0.1 
(solid blue curve), residual oscillations in tc are below lAf, as can be seen from the 
inset of Fig. [31 We therefore recommend a matching width of at least 6uj/lu = 0.1, 
and adopt this value for the remainder of this paper. 

To the lowest PN order, the time to merger is proportional to A/w"*/'^. Therefore, 
a relative frequency width of 10% corresponds to the first quarter of the NR waveform. 
The PN waveform must be sufficiently accurate throughout the entire matching 
interval, so this necessitates a NR waveform which extends considerably into the regime 
where PN errors are small. 

4-3. Systematic errors in the NR waveform 

Let us now examine two errors that arise during a numerical simulation: numerical 
truncation errors and errors introduced by gravitational wave extraction at finite 
radius. 
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We assess the truncation error by constructing hybrids from NR waveforms 
at three different numerical resolutions (N4, N5 and N6 in Rcf. [3^1 [52] with 
approximately 57'^, 62'^, and 67'^ grid points respectively). These hybrids utilize a 
TaylorTS PN waveform and Muj^ = 0.042. We then compute between 
low and medium resolution, and between medium and high resolution. Because of the 
rapid convergence of the spectral methods used to compute these numerical waveforms, 
the errors essentially represent the error of the lower of the two resolutions 

involved. The results are plotted as the dashed-dotted lines in Fig. H] Numerical 
truncation error is seen to be irrelevant, except perhaps at the lowest resolution. This 
result agrees with the results of [13] for the Llama code, and [17] for the BAM code, 
although |27j considers only the numerical waveform and not a hybrid. 

Turning to the gravitational wave extraction at finite radius i?GW: we construct 
hybrid waveforms based on NR waveforms at two finite extraction radii, and one 
where the waveform is extrapolated to infinite extraction radius (we use a TaylorTS 
PN waveform with Mu!,n = 0.042). The solid lines in Fig. H] report the differences 
between the "finite-radius hybrids" and the hybrid based on the extrapolated NR 
waveform. The effects of finite-radius wave-extraction are significant, making even the 
i?GW = 385M numerical simulation unusable for parameter estimation at moderately 
large pcs- Thus, it is imperative to always utilize waveforms extrapolated to infinite 
extraction radius, or to use Cauchy characteristic extraction |57| . This result is of 
relevance to other studies such as [13] , because most of the NR waveforms used in that 
analysis have a finite extraction radius. One should note that the results presented 
here are specific to the numerical waveforms we study, and may not hold in general. 
For example, the results from [SS] show a lower error for finite extraction radius 
waveforms, but were obtained with a different evolution system and gauge (moving 
puncture BSSN, rather than Generalized Harmonic). 

Wave extrapolation to infinite extraction radius is already well established [^[55] . 
so one might view an investigation of the impact of finite Rqw as less important. 
However, we can use the finite Rqw waveforms to comment on several topics of 
relevance. 

We begin by investigating the relative importance of phase- and amplitude-errors. 
Using the Fourier-transforms Sh{f) and h{f), we follow Lindblom et al. [52] and 
write Sh{f) = [5x{f) + ^^<P{f)]h{f), so that 6x{f) and S(j){f) represent the fractional 
amplitude-error and the phase-error of the waveforms in the Fourier-domain. Defining 
noise-weighted averages, 6x = \\Sxh\\/\\h\\, Scf) = (where ||(5x/i|| is 

obtained by substituting 5x{f )h{f ) into the integral in Eq. [5]), Lindblom et. al. show 



Fig. [5] plots 5x and 5(j) for three representative comparions considered in this paper. 

In the case of numerical errors the phase errors dominate over the amplitude 
errors, as can be seen in the left panel of Fig. |5J The right panel of Fig. H] shows 
the amplitude and phase errors for a change in PN approximant and for different 
matching frequencies, which we will discuss in more detail in Sec. 14.41 In these 
cases, amplitude errors sometimes are similar to phase errors, especially for very high 
matching frequencies or very large total mass. In general, however, it is fair to make 
the assumption that phase errors dominate over amplitude errors. The discontinuity 
in 5x and 5(j) for one of the plotted camparisons arises because the best-fit alignment 



that 




(28) 
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Figure 5. Contribution of amplitude- and phase— errors to the overall error 
||(5h||/||/i||. Left: the effect of wave-extraction at finite radius Rqw = 385 Af, 
and the effect of numerical truncation error (see Fig. |4||. Right: the effect of a 
change of Post-Newtonian approximant (TaylorTS vs. TayIorT4), and changes to 
the matching region {uim = 0.038 vs 0.049 and 0.064), cf. Sec.BTIl 



between the two hybrids being compared discontinuously jumps between two local 
minima as M passes through a critical value where the global minimum jumps from 
one to the other minimum. 

Next, we recall that numerical simulations arc particularly susceptible to phase 
errors during the inspiral phase, where the energy flux is small and the inspiral 
time-scale is long. Finite-radius gravitational waveforms induce errors with similar 
properties, as can be seen in the left panel of Fig.lH which plots the phase-errors of the 
finite i?GW numerical waveforms. Therefore, finite-extraction waveforms are a good 
model to assess the importance of phase-errors during the inspiral phase. Comparing 
the left panel of Fig. [Hlto Fig. |H we find that phase-errors of ~ 0.5 rad during the 
numerical simulation - whatever their origin - are not indistinguishable in advanced 
detectors and should be avoided. 

The right panels of Fig. IH] demonstrate that the phase-errors and the respective 
values are proportional to l/i?GW- This confirms our earlier assertion 
that ||(5/i||/||/i|l is proportional to the dominant error that enters a hybrid waveform 
(cf. the discussion after Eq. (fTU|)). This proportionality allows us to predict how 
small numerical phase-errors should be, namely < 0.2 rad to be undetectable at 
Peff = 40, and proportionally smaller for higher SNR. To include some safety margin, 
we recommend a target of about 0.1 rad phase error when plotted similarly to Fig. [5] 
as a function of frequency. The tolerable phase error during the inspiral will naturally 
depend on many factors, most notably the length of the NR waveform (i.e., the 
hybridization frequency uJm)- We have not investigated scaling with length of the 
NR waveform, so our recommendation only holds formally for Mu>m ~ 0.04. 

Fig. [7] uses the comparison between finite i?GW and extrapolated waveforms 
to motivate statements that we made in the context of Eq. ((T7)) . We plot M and 
(ll'^^ll/ll^ll) /2- These quantities should agree if /iii is negligible, which indeed is the 
case. 

Finally, Fig. [7] also shows mismatches between the pure NR waveforms (not 
hybridized) at finite and infinite extraction radius. These mismatches are only 
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Figure 6. Left: Phase difference between the NR waveforms extracted at finite 
-Rgw relative to a waveform extrapolated to infinite extraction radius. Right: 
"Convergence test" of finite Rgw waveforms with extraction radius. The top 
right panel shows the phase error, the bottom panel and both are seen 

to decay like l/if,Qw 



meaningful for sufficiently high mass M so that the NR waveform begins at frequencies, 
/iNR, below the Advanced LIGO frequency band at lOHz, i.e. for M > 130Mq. As 
expected, we see that for M > ISOMq the pure NR mismatches agree very well 
with the hybrid mismatches, because hybridization is unnecessary and unimportant 
for such high masses. However, at lower masses the hybrids have much larger errors 
than at higher masses. The mismatch, for instance, reaches a maximum 10 times 
larger than the maximum mismatch of the NR-only waveforms. Therefore, one can 
use pure NR mismatches only to ascertain the usability of NR waveforms for high 
masses. Small mismatches of pure NR waveforms at high masses do not imply small 
mismatches of hybrid waveforms at lower masses. The only way to assess the quality 
of NR waveforms that are intended for use in hybrids lies in the construction of hybrids 
and computing errors based on these hybrids. Pure NR mismatches are commonly 
used because of their convenience (e.g., in the Samurai project but one must 

keep their limitations in mind. 

All hybrids used in this section use the TaylorTS PN approximant. We have 
repeated the calculations with other approximants without any change in our results. 
Our findings are, therefore, entirely due to the behavior of the numerical waveforms. 

4- 4- Choice of matching frequency 

This section addresses at what frequency hybridization should occur, i.e., how long 
the NR simulations must be. There is a trade-off between computational expense and 
accuracy. On the one hand, we wish to match PN and NR waveforms as close to 
merger as possible to reduce the length and computational cost of NR simulations. 
On the other hand, PN waveforms diminish in accuracy toward merger. Because of 
its importance, this is a very active research topic [M j [27 l [28 ] . 

As a first step towards understanding the importance of the matching frequency 
ojm , we consider a series of hybrids where TaylorTS is matched to the NR simulation 
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Figure 7. Verification of Eq. I|15p relating and mismatch. For 

the finite radius hybrids shown in Figs. [4] and |6l we plot mismatch M and 
(||(5h||/||/i||)^/2, which are virtually identical. We also plot mismatches between 
pure NR waveforms (no hybrids) for so high masses that the initial GW frequency 
of the NR waveform is /it-jpi > lOHz. 

at different values uJm, using a matching interval 6u! = O.lWm- We then compute 
the error measure between a "reference hybrid" which is matched at the 

lowest possible frequency {Mcom = 0.038), and "trial hybrids", matched at higher 
frequencies. The results are presented in the left panel of Fig. [H With increasing 
ojm of the trial hybrid, ||(5/i|j/||/i|| becomes unacceptably large. With decreasing w™, 

decreases, as one would expect, and for Mwtriai = 0.042, the error criterion 
appears to be satisfied even at poff = 40. However, this decrease in ||(5/i||/||/i|| is a 
combination of two different origins. First, with earlier matching, the trial hybrid 
will be closer to the exact waveform, obtainable if one could match arbitrarily early, 
— ^ 0. This is indeed the effect that we arc trying to measure. Second, as the 
matching frequencies of trial- and reference-hybrid approach each other, the hybrid- 
waveforms themselves also approach each other, and |j(5/i||/|i/i|| will decrease. To allow 
conclusions about the trial hybrid, we must quantify the importance of this artificial 
suppression of ||(5/i||/||/i||. 

The right panel of Fig. [8] attempts to determine the importance of the artificial 
suppression due to close matching frequencies. We fix wtriai, and consider decreasing 
values of Wrof- In the limit coj-at — >■ 0, the difference ||(5/i||/||/i|| will measure the error 
in the trial hybrid. To reduce the dimensionality of the plot, the right panel of Fig. [5] 
shows the maximum of over the mass-range [5Mq,200Mq]. As expected, 

becomes smaller and approaches as Wrcf — >■ Wtriai- As was just alluded to, 
this is the artificial suppression of Let us now consider the opposite limit of 

Wref 0. For large wtriai (for instance the top-most, magenta curve {MtJtnai = 0.064)), 
||(5/i||/||/i|j indeed levels off as uti-ci is made as small as possible, and asymptotes 
towards what appears to be some maximum value. We interpret this maximum 
value to represent the true error in the trial hybrid waveform if our NR waveform 
were arbitrarily long. For A/wtriai = 0.064, this asymptotic value is « 0.5, 

demonstrating that a TaylorTS hyrid matched at Mujm — A/wtriai = 0.064 is not 
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Figure 8. Errors of TaylorTS hybrids, as a function of matcliing frequency uJrn- 
Left: ||<5h.||/||/i|| between matching frequencies MiJtrial and fixed Muij.^{ = 0.038. 
The dotted lines indicate upper hmit of ||(5h||/|jh|| for two different effective signal 
to noise ratios p^ff. Right: ||(5/i||/||/i|| maximized over M, as a function of the 
matching frequency of the reference hybrid for trial hybrids using several values 
of (Atrial- The dashed lines are the maximum values of ||(5/i||/||/i|| for each trial 
hybrid compared to a TaylorT4-hybrid matched at Mujm = 0.042. 

useful for any data-analysis purposes. 

As we move toward smaller Wtriai, the error falls. Unfortunately, once 

Atrial is so Small that the error begins to get interesting (e.g. the lowest 

two curves, Mwtriai = 0.046,0.049), the asymptotic behavior as Wrcf —J- is no longer 
visible. In order to re-establish this behavior, we would require such a small value 
in oJrcf such that it is below the starting frequency of the numerical waveform. By 
utilizing only a single PN approximant, we can, therefore, learn only the following: the 
15-orbit SpEC waveform is long enough to ascertain that matching at high frequencies 
Muj„i > 0.055 results in hybrids with unacceptably large errors. However, the NR 
waveform is too short to establish reliably the error incurred by matching at lower 
frequencies 0.038 < A/w,„ < 0.055. 

To gain more insight into the effects shown in Fig.|8l we explore hybrid waveforms 
using different PN Taylor approximants. A natural starting point involves constructing 
our reference hybrid using a PN TaylorT4 approximant, which is known to be 
considerably more accurate than TaylorTS in the non-spinning equal mass case (see 
Ref. [H] for detailed discussion) . If we assume that TaylorT4 incurs a negligible error 
relative to the "true" inspiral waveform, error criteria between TaylorT4 

hybrids and the previously considered TaylorTS hybrids will display the error in the 
TaylorTS hybrid. These errors are indicated by the horizontal dashed lines in the 
right panel of Fig. HI (Because TaylorT4 agrees so well with NR, for convenience 
here we use only one TaylorT4 hybrid matched at Mlj = 0.042). These dashed 
horizontal lines in the right panel of Fig. \E\ are consistent with our conclusions of the 
preceeding paragraph: the errors in the TaylorTS hybrids matched at high frequency 
(Mwtriai = 0.058, 0.064) can be resolved by pushing uji-d to the start of the NR 
waveform, and the asymptotic value of the TaylorTS-only comparisons agrees with 
those of TaylorTS — TaylorT4. For smaller Wtriai, we cannot reach the asymptotic 
regime, and while ||(5/7,||/||/ij| for TaylorTS-only comparisons increases as Wrof is pushed 
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Figure 9. Left: Errors of TaylorTl hybrids as a function of matching frequency. 
Right: Errors of TaylorT4 hybrids as a function of matching frequency. These 
plots mirror the right panel of Fig. |8l except using a different PN approximant 
(note the different y-axis scale in the right panel). 



as low as possible, the error is still far away from the expected limit given by the 
TaylorT3-TaylorT4 comparison. The dashed TaylorT3 — TaylorT4 lines indicate errors 
too large even for peff = 10. Under our conservative error criteria, we conclude that 
the NR waveform is too short to construct useful TaylorTS hybrids. 

So far, we have used the TaylorTS PN waveform, because its analytical nature 
makes it the simplest to implement. We now consider different PN approximants. 
The left panel of Fig. [H plots the same information as the right panel of Fig. [HI 
except that TaylorT3 is now replaced by TaylorTl. Similar to the right panel of 
Fig. [3 we can clearly see the suppression of the differences when Wtriai is too close to 
Wrof- As above, this effect manifests itself in oc jwrof — i^triail- However, as 

opposed to TaylorT3, for the TaylorTl approximants, we do not see any indication 
that ||i5/i||/||/i|| levels off as Wrof is made smaller, not even for the largest possible 
difference (wtriai ~ Wrof )• Our conclusions for TaylorTl hybrids, therefore, are identical 
to those reached for TaylorT3 hybrids; the current 15-orbit NR waveform is too 
short to construct reliable TaylorTl hybrids on the basis of our error criteria. Or 
conversely, TaylorTl is too inaccurate to be utilized with the currently available 
numerical waveforms within our error framework. 

The right panel of Fig. H] shows the same comparisons using the TaylorT4 PN 
approximant. TaylorT4 is special, in that it happens to agree very well with the 
numerical simulations for the equal-mass, non-spinning case. Indeed, the right panel 
of Fig. ini confirms this fact; the errors arc far smaller, within our limits 

for detection and parameter estimation even for poff = 100 for matching frequencies 
Mujm < 0.052. For all Wtriai, convergence to an asymptotic value is visible as Wref is 
decreased. Therefore, TaylorT4 hybrids are suitable in the equal-mass, non-spinning 
case for event detection and parameter estimation. However, as is well known, this 
agreement is coincidental and does not carry over to more generic configurations |60) ; 
the exceptional good behavior of TaylorT4 is, thus, not helpful for the generic case. 
Nevertheless, convergence towards an upper error limit at lower matching frequency 
suggests the suitability of different PN approximants at frequencies lower than current 
NR waveforms can provide. 
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Figure 10. Error criterion evaluated between TaylorTS hybrids at different post- 
Newtonian order. All hybrids matched at Muim = 0.042, and all comparison 
relative to a hybrid using Taylor T3 3.5/3.0 (3.5PN order in phase, 3.0 PN order 
in amplitude). 



We point out that our results presented in Figs. [S] and [5] are consistent with the 
results of Buonnano et al. [M| who compute overlaps for M = 20Mq between non- 
hybridized Taylor approximants terminated at ISCO. Translating into our notation, 
they find ||(5/i||/||/i|| = 0.64 when comparing Taylor T4 to Taylor T3, and ||(5/i|j/||/i|| = 
0.34 when comparing Taylor T4 to Taylor Tl. At our highest matching frequency 
Muitriai = 0.064 (which is closest to pure PN waveforms), one can read off from 
Figs.ilandinithat = 0.56 and |l(5/i||/||/i|| = 0.46 respectively. This agreement 

is reasonably good, considering the differences in the examined waveforms. 

To close this section, we briefly investigate hybrids constructed from TaylorTS at 
different PN orders. We construct all hybrids at a matching frequency of MtOm = 0.042 
and compute between a reference hybrid using the most accurate PN order 

(3.5PN in phase, 3PN in amplitude) and trial hybrids with lower PN expansion orders. 
The legend of Fig. (TU] indicates which PN orders we consider. All differences, even 
between 3.0 and 3.5PN phase accuracy, are so large, suggesting that lower order PN 
approximants may not be sufficiently accurate for data-analysis purposes. This finding 
is consistent with our earlier findings; we established earlier that using our error 
criteria, TaylorT3 3.5/3.0 is insufficient when matched at Moj„i = 0.042. As Fig. [TO] 
confirms, lower-order PN approximants are even less accurate, and consequently, are 
also insufficient within our analysis. 

4-5. Effects of the noise curve 

We have so far considered detection and measurement accuracies only in the context 
of a single Advanced LIGO interferometer. This section briefly presents ||(5/i||/||/ij 
computations for different instrument noise curves including Initial LIGO (with a low- 
frequency cut-off of 40IIz) and a frequently used analytical fit [13] to the Advanced 
LIGO sensitivity curves. These noise curves are plotted in Fig. [T] The overall scale 
of the noise curve Sn{f) and the overall scale of the waveforms h{t) cancels in the 
normalized error criterion Nevertheless, the shape of the waveforms and of 

Sn{f) are important |28| . Specifically, a noise-spectrum with a wider bandwidth will 
be sensitive to a larger portion of the waveform. 
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Figure 11. Effect of different noise curves on ||<5h||/||/i||. Shown are differences 
between finite extraction radius i?GW = 385A/ and waveforms extrapolated to 
infinity, cf. Fig. |4] 



In order to investigate the impact of different noise curves, we use the waveform at 
finite extraction radius i?GW = 385A/, already presented in Figs.UHZl Fig. [TT] repeats 
this earher calculation for all noise curves presented in Fig. [U The error 
is larger for the larger-bandwidth noise curves, in particular ZERO_DET_HIGH_P and 
ZERD_DET_LOW_P. Comparing ZERO_DETXOW_P and ZERD_DET_HIGH_P, we note that 
ZERD_DET_LDW_P is more sensitive at low frequencies relative to high frequencies, and 
correspondingly it results in a larger at high masses (i.e., low- frequency 

waveforms) and smaller at low masses, compared to ZERO_DET_HIGH_P. 

We also include the ZERO_DET_HIGH_P noise curve with a low-frequency cutoff at 20 
Hz. This increased low-frequency cutoff has no impact on low mass systems with 
M < 2OM0, but results in reduced sensitivity (i.e. smaller ||(5ft.||/||/i||) at larger masses, 
compared to the unmodified ZERO_DET_HIGH_P noise curve. 

We note that the analytical fit to the Advanced LIGO noise curve (green line) 
underestimates the error compared to the projected Advanced LIGO noise 

curve. Also, hybrid waveforms, sufficient for measurement accuracy with Initial LIGO, 
may no longer meet the same measurement criteria within an identical mass range and 
Peff in the case of Advanced LIGO, as seen in Ref. [3S]. Besides the larger bandwidth 
(which directly affects ||5/i||/||/i||), the Advanced LIGO noise curves are more sensitive 
than Initial LIGO. This needs to be taken into account when choosing an acceptable 
PeS- The effect of noise curve has also been studied in Ref. [M], where the error 
between different resolutions of NR waveforms and different types of NR codes is 
calculated through the Advanced LIGO and Initial LIGO noise curves. A similar 
trend is observed with respect to an increase in sensitivity (and therefore error) at 
high masses. 

5. Discussion 

In this paper wc have presented a comprehensive study of errors that affect hybrid 
waveforms. We assess the quality of the hybrids based on their suitability for Advanced 
LIGO parameter estimation, using the criterion most recently presented in [52], and 
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refined in |28| . 

ll^ll Poff 

The left-hand side of this inequality is independent of the source distance, and 
therefore, forms a convenient quantity for this analysis. The right-hand side 
incorporates the single detector signal-to-noise ratio p of to-be-analyzed events, a 
safety factor e [28] and possibly a correction factor to account for a network of 
detectors, cf. Eq. PT|) . The inequality provides a bound on p (or, cquivalcntly, 
distance) below which the error 5h is undetectable in the GW detectors, indication 
that h is sufhcicntly accurate. As argued in Sec. [51 poff = 40 is a reasonable value, 
on which the following conclusions are based. However, one should keep in mind that 
Peff = 100 — 250 might be relevant if optimistic assumptions about the number of 
GW sources turn out to be correct, resulting in stronger accuracy demands on hybrid 
waveforms. 

We investigated a large variety of effects that cause possible errors 5h. The 
general theme of our findings is that Eq. ([2^ places strong accuracy requirements on 
all elements of the hybrid construction: 

(i) Section B?T] shows that the PN waveform has to be aligned correctly to 5tc < IM 
time-offset, about 1/100 of a gravitational wave cycle. 

(ii) Section 14.31 shows that the numerical waveform must have phase-errors < 0.1 
radians. 

(iii) Section 14.41 indicates that PN Taylor-approximants in the time domain are not 
sufficiently accurate when representing the inspiral waveform up to the start of 
currently available numerical relativity waveforms. 



These three findings are interrelated. We have already pointed out in Eq. ([27]) 
that the 3.5PN order contributes a time-offset of about <3.5pn ~ lOOM at typical 
matching frequencies Muj k, 0.04. The analytical formulae from the TaylorT2 
approximant (211 [U] show that the 3.5PN order for an equal-mass binary contributes 

03.5PN = -\^^^i^{Amf'^ « -2.6 rad (30) 

to the GW phase (at orbital frequency MVl Mbj/2 = 0.02). Both the PN 
contributions ^s.spn and (/js.spn arc much larger than the limits on the time and phase 
errors we established for the numerical waveforms (lAf and 0.1 rad, respectively), 
so it is unsurprising that the time-domain Taylor approximants are insufficient, 
or conversely, that the numerical waveforms are too short within our strict error 
requirements. 

So how long do NR waveforms need to be in order to reliably attach a Taylor-PN 
inspiral? Or equivalently, how early must one stop using the PN waveforms? Fig. [T2l 
suggests a conservative answer to this question. Both panels of Fig. [T^j show the same 
data, but with different axes. The data presented arc for M = lOM© and 

M = 20Afo between pairs of {TaylorTl, TaylorT3, TaylorT4} hybrids, aU matched 
at the same frequency Wm, and plotted as a function of hybridization frequency a;,„. 
Because the TaylorT hybrids only differ in their treatment of uncontrolled higher-order 
PN terms, these differences are a reasonable measure of the truncation error of the 
PN series. Fig. [T^] shows that these differences are very large, far above the limit of 
Eq. ([2^. The right panel of Fig. [T^ attempts power-law fits to as a function 
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Figure 12. Differences between hybrids using different PN approximants, but 
matched at the same matching frequency. Plotted is the value of ||(5/i||/||fc|| at 
IOMq and 20Mq, as a function of u)rn- Both plots show the same data, albeit 
with different axes. The right plot contains power law fits to the smallest four 
data-points. 



of LOm- Extrapolation based on fits is liighly unreliable, but nevertheless, this figure 
indicates that one must match at significantly lower frequencies Wm than currently 
available. The time-to- merger scales in proportion to (A/w)"^/"^. so this implies a 
need for NR waveforms several times longer than those currently computed. 

Figure [12] shows data for total mass of \QMq and 2QMq, with the lower total mass 
resulting in larger Our findings are compatible with the results of earlier 

PN-inspiral and EOB work [55]. In the case of the significantly lower mass (1.42, 
1.38) Mq, using their faithfulness quantity, Buonnano et al. show that Taylor T2, 
T3, and T4 exhibit a corresponding j|(5/i||/||/i|| value of 0.14 with each other, whereas 
Taylor Tl has a somewhat worse agreement of 0.24. However, we reiterate that the 
criterion p9|) is sufficient, and may therefore place unnecessary strong constraints on 
waveform accuracy. 

Figure [T^ clear Iv presents a challenge for the current state of the art in numerical 
and analytical waveform modelling. Fortunately, there are several potential avenues to 
reduce the "frequency gap" (in the language of ^E\) between PN and NR waverforms. 

First, one can compute longer NR waveforms. This is clearly challenging, as the 
following considerations based on the SpEC code illustrate. By order of magnitude, a 
simulation of a generic inspiral with moderate mass-ratio and spins currently requires 
several 10^ CPU-hours, lasting several months. These numbers apply to an evolution 
time of about 5000A/ at accuracies indicated by our analysis here. Running SpEC for 
~ 30 orbits is feasible, however, this will require wall-clock times approaching a year. 
We caution that no careful convergence tests have been performed for such simulations, 
though preliminary simulations appear promising. Excessive wall-clock time is the 
most restrictive factor in simulations longer than 15 orbits. This situation can 
be ameliorated with at least three different approaches: increase parallelism to many 
more cores; utilize a faster computing architecture, e.g. graphics accelerators |61j : 
or develop more efficient numerical algorithms (HU [32 [33]. All three approaches 
are promising, and we expect that a combination of them will result in a significant 
extension of computational capabilities. 

An alternative approach is to improve the accuracy of the PN series, either 
through calculation of higher PN order, or by seeking more rapidly converging 
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analytical representations of the inspiral phase. The EOB formalism provides such 
a representation, as recently argued in |28j . However, some improvements to the 
EOB formalism (e.g., [62]) were developed after high-accuracy NR waveforms became 
available for comparisons. While the improvements are well motivated on physical 
grounds, the additional confidence that EOB models are developed independently from 
the numerical results no longer applies. Finally, once one considers tuned models, one 
loses all ability to use differences (e.g., to the numerical result) as an error measure. 
Such models, by definition, will agree very well in the frequency range in which they 
are tuned to numerical simulations. However, such tuning does not guarantee that the 
cycles before the tuning region are well represented by the tuned model. Even with 
these reservations, EOB models are perhaps the most promising avenue to close the 
"frequency gap" between numerical simulations and Taylor post-Newtonian models. 

Finally, a more detailed analysis of waveform accuracies for GW data- analysis 
might result in weakened requirements. The impact of detector calibration errors 
or of non-Gaussian detector-noise might overwhelm the limits presented here that 
are based on the ideal detector case. Furthermore, parameter estimation is mostly 
sensitive to waveform errors 6h tangential to the signal manifold. If the errors in 
hybrid-waveforms are predominantly orthogonal to the signal manifold, then Eq. ([2^ 
may be overly restrictive. 

In addition, let us briefly summarize how our findings impact current efforts 
in GW detection. The Ninja-2 project [13 [IB] is well underway. This project has 
chosen fairly loose accuracy requirements to encourage wide participation, namely 
NR waveforms that must be hybridized at a frequency MtUm, < 0.075, with phase 
accuracy of the numerical waveform better than 0.5 rad. This phase accuracy seems 
sufficient for event detection purposes (cf. Figs. ID and [6]), but matching at such high 
frequencies will be dominated by PN errors to a degree that would likely impact 
parameter estimation, (cf. Fig. [T^ . While the Ninja-2 project focuses primarily on 
event detection, some of its efforts are geared at parameter estimation. These efforts 
will be most useful in shedding light on the accuracy requirements discussed in this 
paper. Future Ninja projects may have to sharpen their accuracy requirements to 
avoid biasing the results by inadequate model waveforms and will likely be targeted 
to improve understanding of parameter estimation. 

This work can of course be extended in several directions. Most immediately, 
our conclusions are limited by the length of the employed NR waveform, despite it 
being the longest published NR simulation. Therefore, it will be useful to revisit 
the equal-mass, non-spinning case when longer numerical data become available. 
Furthermore, the study needs to be extended to BBHs with different mass ratios and 
spins. With increasing mass ratio, the binary spends more orbits in the strong-field 
regime (inversely proportional to the symmetric mass ratio) , so it is likely that unequal 
mass binaries will require yet longer numerical simulations than the equal-mass case 
considered here, as also pointed out by [28]. Spin effects are known only to lower PN 
order, so it is likely that black hole spin may necessitate longer NR simulations as 
well, dependent on their initial spin orientation. 

Besides these obvious extensions based on refined numerical simulations, there 
are also conceptual questions that deserve further attention. Eq. ([^^ incorporates 
a detector network only in the crudest possible way. It would be interesting to 
investigate in more detail how Eq. ([29]) generalizes to a network of multiple detectors. 
More detailed investigations into accuracy standards would also be valuable. On the 
one hand, Eq. (|29p is a sufficient condition and necessary accuracy standards for 
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parameter estimation may be considerably weaker. On the other hand, even when 
a certain error Sh satisfies Eq. ([2^]) . it might nevertheless lead to systematic bias in 
parameter estimation [63]. One must also consider how the presence of calibration 
errors affects the present conclusions. Finally, throughout this paper we have applied a 
time and phase shift to minimize \ \6h\ \ between these two waveforms. This is generally 
acceptable, but doing so discards information about the time of merger and phase 
of merger of the binary. A separate study will be necessary to determine how to 
measure the binary's time of merger as accurately as possible, for instance, in studies 
of electromagnetic counterparts. 
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